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INTRODUCTION 


Breast  cancer  has  been  considered  as  a  major  problem  and  the  most  common  cancer  among 
women.  Annually,  a  total  of  348,000  cases  of  breast  cancer  is  diagnosed  and  almost  1 15,000  are 
killed  by  this  in  the  US  and  European  Community  [1l  Tremendous  efforts  have  been  made  in  the 
incremental  improvements  in  imaging  technologies  in  the  field  of  breast  cancer  detection. 
Mammography  is  currently  the  most  important  and  efficacious  tool  for  the  early  detection  of 
breast  cancer  [2].  However,  the  nature  of  the  two-dimensional  mammography  makes  it  very 
difficult  to  distinguish  a  cancer  from  overlying  breast  tissues,  especially  for  dense  breast  cases. 

Digital  breast  tomosynthesis  is  a  three-dimensional  breast  imaging  method  that  allows  the 
reconstruction  of  an  arbitrary  set  of  planes  in  the  breast  from  limited-angle  series  of  projection 
images  as  the  x-ray  source  moves  along  an  arc  above  the  breast.  A  variety  of  tomosynthesis 
reconstruction  algorithms  have  been  proposed  including  the  traditional  shift-and-add  (SAA),  the 
image-stretching  method  proposed  by  Nilklason  and  colleagues  the  maximum  likelihood 
iterative  algorithm  (MLEM)  by  Wu  et  al.  [4’5-',  tuned-aperture  computed  tomography  (TACT) 
reconstruction  methods  developed  by  Webber  and  investigated  by  Suryanarayanan  et  al  [6,7], 
algebraic  reconstruction  techniques  (ART)  [8'9’10-',  filtered  back  projection  (FBP)  I11-12’13],  an(j 
matrix  inversion  tomosynthesis  (MITS)  [14'15]. 

The  purpose  of  this  project  is  to  optimize  and  compare  several  different  tomosynthesis 
reconstruction  algorithms  that  are  either  initially  investigated  in  our  lab  or  currently  very 
popular,  and  to  optimize  the  imaging  acquisition  techniques.  Based  on  this  project,  we  hope  to 
contribute  to  the  optimal  breast  tomosynthesis  technique  for  better  breast  cancer  detection. 

BODY 

Task  1.  Optimization  of  different  candidate  tomosynthesis  reconstruction  methods  (Month 
1-12): 

1,1.  Select  candidate  algorithms  from  different  algorithm  families  for  optimization  and 
comparison.  2-4  candidate  algorithms  from  different  algorithm  families  will  be  chosen. 
(Month  1-4). 

We  have  investigated  a  few  reconstruction  algorithms,  including  shift-and-add  (SAA)  algorithm, 
Nilklason’ s  image  stretching  shift-and-add  (NIKL),  maximum  likelihood  expectation 
maximization  (MLEM),  matrix  inversion  tomosynthesis  (MITS),  and  filtered  back  projection 
(FBP). 

We  found  that  the  traditional  shift-and-add  (SAA)  is  a  common  mathematical  method  to  line  up 
each  projection  image  based  on  its  shifting  amount  to  generate  reconstruction  slices.  MITS  was 
originally  invented  in  our  lab  by  Dobbins  [15]  and  we  applied  it  to  the  breast  tomosynthesis 
imaging  successfully  [14l  MITS  shows  better  high  frequency  response  in  removing  out-of-plane 
blur.  We  also  developed  our  FBP  reconstruction  based  on  central  slice  theorem  and  Fourier 
frequency  sampling  density.  In  order  to  control  the  high  frequency  noise  amplification, 
Hamming  Gaussian  filters  were  designed  and  applied  to  our  FBP  algorithm  ri2l  MLEM 
algorithm  is  an  effective  iterative  method  in  breast  tomosynthesis  reconstruction.  However,  it  is 
time-consuming  due  to  intensive  computation.  Therefore,  in  this  project,  we  selected  SAA, 
MITS  and  FBP  as  our  three  candidate  algorithms  for  comparison  and  optimization.  But  we  also 
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did  related  research  on  other  algorithms  such  as  point-by-point  back  projection  (BP),  NIKL,  and 
MLEM. 

Work  for  this  specific  task  also  resulted  in  part  of  a  manuscript  that  we  submitted  to  the  journal 
of  Medical  Physics  for  publication  (reportable  outcome  #1).  During  our  investigation  of  different 
algorithms,  we  found  that  quite  a  few  other  algorithms  depend  on  a  traditional  SAA  method. 
However,  traditional  SAA  is  not  appropriate  for  the  isocentric  motion  in  digital  breast 
tomosynthesis  because  it  doesn’t  take  into  the  account  of  shift  amount  at  the  direction  orthogonal 
to  tube’s  motion  direction.  A  simple  SAA  reconstruction  algorithm  is  not  entirely  suitable  for 
breast  tomosynthesis,  especially  for  small  structures  such  as  microcalcifications,  which  have  an 
important  bearing  in  clinical  breast  cancer  detection  tasks.  The  manuscript  we  submitted  focused 
on  the  importance  of  point-by-point  back  projection  for  reconstruction  of  small  structures  such  as 
microcalcifications. 

1.2.  Characterize  the  effect  of  three  acquisition  parameters  including  total  Tomographic- 
Angle  (TA),  Number  of  projection  images  (N),  and  Reconstruction-Slice-Spacing  (RSS)  for 
each  reconstruction  algorithm,  according  to  physical  measurements  of  impulse  response 
analysis,  modulation  transfer  function  (MTF)  and  noise  power  spectrum  (NPS).  (Months  5- 
12). 

1.2.1.  Simulate  impulses  with  different  acquisition  parameters.  Apply  each  candidate 
tomosynthesis  algorithms  to  generate  slice  reconstruction  images.  Analyze  the  impulse 
response  of  each  candidate  algorithm.  (Months  5-7) 

Work  on  this  task  began  well  ahead  of  schedule  during  the  first  year.  We  simulated  impulses  by 
ray-tracing  method  with  a  few  different  combinations  of  acquisition  parameters.  The  candidate 
tomosynthesis  algorithms  were  used  to  generate  slice  reconstruction  images.  The  impulse 
response  of  each  candidate  algorithm  was  analyzed. 

Parameters  of  a  selenium-based  direct  conversion  Siemens  Mammomat  Novation  DR  prototype 
system  was  modified  to  be  used  for  geometries  of  the  simulation.  Two  different  impulse 
locations  were  investigated  in  the  project:  1)  an  impulse  that  is  exactly  underneath  the  x-ray 
source  (near  the  chest  wall)  and  in  a  defined  reconstructed  plane  (20  mm  above  the  detector 
surface);  2)  an  impulse  that  is  approximately  4  cm  away  from  the  chest  wall  and  in  a  defined 
reconstructed  plane  (20  mm  above  the  detector  surface).  Datasets  of  the  impulse  with  different 
number  of  projections  images  N  and  total  Tomographic  Angle  TA  of  the  simulated  x-ray  point 
source  were  simulated  and  reconstructed  by  each  candidate  algorithm  for  comparison. 

We  found  that  the  shift-and-add  method  (SAA)  was  similar  to  the  Niklason  method  (NIKL)  in  all 
cases  examined.  Compared  with  SAA  and  NIKL  methods,  MITS  and  FBP  were  always  better  for 
removal  of  out-of-plane  blur  artifacts.  For  optimization  of  Number  of  projection  images  (N),  not 
much  difference  was  noted  between  cases  with  different  N  for  comparison  of  SAA.  For  MITS 
and  FBP,  bigger  N  cases  were  better  in  high  frequency  noise  removal. 

Figure  1  shows  an  example  of  our  investigation  to  optimize  the  N  parameter  when  the  impulse 
was  located  exactly  underneath  the  x-ray  tube  and  halfway  between  two  neighboring  planes, 
which  are  represented  by  Z-0.5mm  and  Z+0.5mm.  In  this  simulation,  Z=19.5mm  above  the 
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detector.  The  Z- 1.5mm  and  Z+1.5mm  locations  are  planes  that  are  1.5mm  lower  and  1.5mm 
higher  than  the  location  of  the  impulse.  In  figure  1,  the  total  tomographic-angle  TA=20°  and  the 
reconstruction-slice-spacing  RSS=lmm.  The  x  axis  represents  the  pixel  location  in  the  column 
containing  the  impulse,  and  the  y  axis  is  the  normalized  impulse  response’s  amplitude.  In  figure 
1,  (a),  (b),  and  (c)  represent  three  selected  candidate  algorithms  of  SAA,  MITS,  and  FBP 
respectively.  We  found  that  there  was  no  substantial  difference  for  11,  21  and  41  projection 
numbers  with  SAA.  With  a  larger  N,  the  MITS  preformed  much  better  and  clearly  reduced  more 
high  frequency  oscillation,  and  provided  clearer  structure.  FBP  preformed  slightly  better  in  out- 
of-plane  structure  removal  at  the  larger  N  due  to  better  sampling  in  frequency  space. 


N=1 1 


(a)  N=21 


N=41 


Z-  1.5mm 


Z-0.5mm 


Z+0.5mm 


Z+ 1.5mm 


N=1 1 


(b)  N=21 


N=41 


Z- 1.5mm 


Z-0.5mm 


Z+0.5mm 


Z+ 1.5mm 


N=1 1 


V  -s-“ 


(c)  N=21 


N=41 


Z- 1.5mm 


Z-0.5mm 


Z+0.5mm 


Z+ 1.5mm 


Figure  1.  Out-of-plane  impulse  response,  near  chest  wall:  (a)  SAA  (b)  MITS  (c)  FBP 
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Work  for  the  specific  task  1.2.1  resulted  in  a  proceedings  abstract  at  AAPM,  a  scientific 
conference  for  medical  imaging  in  2006.  In  the  abstract,  we  used  the  impulse  response  analysis 
method  to  demonstrate  the  improved  reconstruction  of  two-dimensional  SAA  compared  with 
traditional  SAA  algorithm  (reportable  outcome  #2). 

1.2.2.  Analyze  the  MTF  curves  of  different  acquisition  parameters  for  each  candidate 
algorithm.  (Months  8-9) 

During  the  investigation,  we  recognized  that  the  MTF  measurement  actually  includes  two  parts: 
1)  the  system  MTF  of  the  detector;  2)  the  reconstruction  MTF  associated  with  specific 
reconstruction  algorithm  and  acquisition  parameters.  The  system  MTF  describes  the  measured 
MTF  of  the  detector.  A  previously  published  edge  method  L  ’  J  was  applied  at  a  range  of  tube 
angles  to  see  if  there  is  any  difference  with  angle.  Figure  2  shows  the  measured  system  MTF 
when  X-ray  tube  is  located  at  different  angular  locations.  The  MTF  varied  little  with  different 
angles.  The  averaged  MTF  was  calculated  and  served  as  the  system  MTF  of  the  detector. 


The  reconstruction  MTF  was  calculated  as  the  Fourier  transform  of  the  impulse  response  on  the 
reconstruction  plane  associated  with  specific  reconstruction  algorithm  and  acquisition 
parameters.  Figure  3  shows  an  example  how  we  calculated  the  reconstruction  MTF.  A  data  set  of 
tomosynthesis  projection  images  of  a  delta  function  at  40  mm  above  the  detector  surface  plate 
and  40  mm  away  the  chest  wall  was  simulated  with  acquisition  parameters  of  N=25,  TA=25°  and 
RSS=1  mm.  The  simulated  tomosynthesis  sequence  was  reconstructed  by  point-by-point  back 
projection  (BP)  algorithm.  Figure  3(a)  shows  the  impulse  response  where  x  and  y  axes  represents 
the  pixel  location  on  the  image.  Figure  3(b)  represents  the  reconstruction  MTF.  Figure  4  shows 
the  total  MTF  as  the  combination  of  the  measured  system  MTF  and  calculated  reconstruction 
MTF  along  tube’s  motion  direction  and  the  direction  orthogonal  to  tube’s  motion. 
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Figure  3.  Point-by-point  BP  with  N=25  and  TA=25° :  (a)  impulse  response;  (b)  MTF 


Frequency  (cycles/mm) 


cycles/mm 


(a) 


(b) 


Figure  4.  Reconstruction,  system,  and  total  MTFs:  (a)  direction  orthogonal  to  tube’s  motion  direction  (b)  tube’s  motion  direction 


We  applied  the  impulse  response  and  MTF  analysis  methods  in  section  1.2.1  and  1.2.2  to  a 
comparison  of  BP  and  SAA  algorithms.  It  resulted  in  part  of  a  manuscript  we  submitted  to 
Medical  Physics  for  peer-reviewed  journal  publication  (reportable  outcome  #1). 

1.2.3  Analyze  the  NPS  curves  of  different  acquisition  parameters  for  each  candidate 
algorithm.  (Months  10-12) 

We  did  the  NPS  measurement  and  analysis  of  selected  candidate  algorithms  of  SAA,  MITS  and 
FBP  for  comparison  and  optimization.  We  acquired  tomosynthesis  sequences  of  flat  images  with 
five  different  imaging  acquisition  techniques:  1)  N=13,  TA=50°;  2)  N=13,  TA=25°;  3)  N=25, 
TA=50°;  4)  N=25,  TA=25°;  5)  N=49,  TA=50°.  The  total  exposures  were  same  for  tomosynthesis 
sequences  with  different  acquisition  techniques.  Three  reconstruction-slice-spacings  (RSS)  of 
1mm,  2mm,  and  4mm  were  used  for  comparison.  We  also  simulated  and  reconstructed  projection 
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images  of  a  two-dimensional  square  to  match  gains  and  offsets  for  different  candidate 
reconstruction  algorithms  and  imaging  technique  parameters.  Our  previously  published  methods 
[19]  were  applied  to  calculate  the  NPS.  Figure  5  shows  an  example  of  our  NPS  analysis  for  MITS 
algorithm  with  different  acquisition  techniques  with  RSS=lmm. 


Figure  5:  MITS  reconstructions  with  1mm  RSS:  Comparison  for  different  acquisition 
parameters  of  Number  of  projections  (N)  and  Tomographic -Angle  (TA) 


We  found  that  with  the  same  acquisition  parameters  of  N=49  and  TA=50°,  the  SAA  and  FBP 
showed  no  difference  when  we  varied  the  RSS.  The  NPS  from  FBP  reconstructions  showed 
reduced  NPS  at  high  frequencies  due  to  the  applied  low-pass  Hamming  and  Gaussian  filters.  For 
MITS  reconstruction,  compared  with  2mm  and  4mm  slice  spacing,  lnnn  slice  spacing  has  better 
noise  response  at  low-to-middle  frequency  range.  At  middle  to  high  frequency  range,  their  noise 
responses  were  similar.  These  are  also  true  for  other  four  acquisition  parameters  of  N=25  and 
TA=50°,  N=25  and  TA=25°,  N=13  and  TA=50°,  N=13  and  TA=25°. 

Work  for  the  specific  task  1.2  resulted  in  a  proceedings  paper  at  SPIE,  the  primary  scientific 
conference  for  medical  imaging  in  2007  (reportable  outcome  #3).  We  realized  that  the  NPS  itself 
is  not  an  appropriate  method  for  comparison  of  different  algorithms  due  to  different  resolution  of 
each  algorithm.  In  the  paper,  we  proposed  a  methodology  of  noise-equivalent  quanta  NEQ  (f) 
analysis  as  a  better  way  to  compare  and  optimize  different  algorithms  and  acquisition 
parameters. 


Task  2.  Comparison  of  different  candidate  tomosynthesis  reconstruction  methods  (Months 
13-24): 

Generally  speaking,  from  impulse  response,  MTF,  and  NPS  analysis,  we  found  that  MITS  and 
FBP  perfonned  better  with  bigger  N  and  wider  TA  (such  as  N=49  and  TA=50°).  MITS 
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performed  better  with  smaller  RSS  (RSS=lmm).  We  also  found  from  experiments  that  N=25, 
TA=50°,  RSS=lmm  is  an  efficient  way  in  clinical  application  for  MITS  and  FBP.  We  have 
proposed  a  methodology  of  NEQ  (f)  analysis  to  compare  and  optimize  reconstruction  algorithms 
and  acquisition  techniques.  We  are  working  on  the  NEQ  (f)  analysis  now  and  expect  to  finish  the 
NEQ  (f)  analysis  in  May  2007.  Once  technique  parameters  have  been  optimized,  we  will  conduct 
the  comparison  of  different  candidate  tomosynthesis  reconstruction  methods  based  on  the  lesion 
simulation  and  human  observer  study. 


KEY  RESEARCH  ACCOMPLISHMENTS 

•  Investigated  several  different  reconstruction  algorithms  for  digital  breast  tomosynthesis, 
including  SAA,  NIKL,  BP,  MITS,  FBP,  MLEM.  Selected  candidate  algorithms  and 
compared  perfonnance  against  each  candidate  algorithms. 

•  Analyzed  impulse  response,  MTF,  and  NPS  by  simulation,  experiments  and 
characterization  to  compare  and  optimize  the  imaging  acquisition  parameters  including 
total  Tomographic- Angle  (TA),  Number  of  projection  images  (N),  and  Reconstruction- 
Slice-Spacing  (RSS). 

•  Investigated  the  importance  of  point-by-point  BP  for  isocentric  motion  in  digital  breast 
tomosynthesis,  especially  for  reconstruction  of  small  structures  such  as 
microcalcifications.  Proposed  a  methodology  of  NEQ  (f)  analysis  for  comparison  and 
optimization. 

REPORTABLE  OUTCOMES 

The  following  manuscripts  and  abstract  are  attached  at  appendices  1 ,  2  and  3  with  the  same 
numbers.  The  names  of  the  fellow  (Chen)  and  mentor  (Dobbins)  are  boldfaced  for  emphasis. 

1.  Y.  Chen,  J.  Y.  Lo,  and  J.  T.  Dobbins  III.,  “Importance  of  point-by-point  Back 
Projection  (BP)  correction  for  isocentric  motion  in  digital  breast  tomosynthesis:  relevance 
to  morphology  of  microcalcifications,”  Med.  Phys,  in  review,  2007. 

2.  Y.  Chen,  J.  Y.  Lo,  J.  T.  Dobbins  III,  “Two-dimensional  Shift-And-Add  Algorithm  for 
Digital  Breast  Tomosynthesis  Reconstruction,”  Med.  Phys.  33  (6),  2001  (2006). 

3.  Y.  Chen,  J.  Y.  Lo,  N.  T.  Ranger,  E.  Samei,  J.  T.  Dobbins  III,  “Methodology  of  NEQ(f) 

analysis  for  optimization  and  comparison  of  digital  breast  tomosynthesis  acquisition 

techniques  and  reconstruction  algorithms,”  Proc.  SPIE  6510,  65101-1,  (2007). 

CONCLUSIONS 

We  have  investigated  different  reconstruction  algorithms  for  digital  breast  tomosynthesis, 
compared  and  optimized  candidate  algorithms.  We  have  analyzed  the  impulse  response,  MTF, 
and  NPS  by  simulation,  experiments,  and  computation  to  compare  and  optimize  candidate 
algorithms  and  different  acquisition  parameters.  We  also  proposed  a  NEQ(f)  analysis  to  better 
compare  and  optimize  reconstruction  algorithms  and  acquisition  techniques.  Future  work  will  be 
done  to  further  compare  the  algorithms  with  optimized  acquisition  parameters  by  lesion 
simulation  and  human  observer  study. 
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ABSTRACT 

40  Digital  breast  tomosynthesis  is  a  three-dimensional  imaging  technique  that  provides  an 
arbitrary  set  of  reconstruction  planes  in  the  breast  from  a  limited-angle  series  of 
projection  images  acquired  while  the  x-ray  tube  moves.  Traditional  Shift- And- Add 
(SAA)  tomosynthesis  reconstruction  is  a  common  mathematical  method  to  line  up  each 
projection  image  based  on  its  shifting  amount  to  generate  reconstruction  slices.  With 
45  parallel-path  geometry  of  tube  motion,  the  path  of  the  tube  lies  in  a  plane  parallel  to  the 
plane  of  the  detector.  The  traditional  SAA  algorithm  gives  shift  amounts  for  each 
projection  image  calculated  only  along  the  direction  of  x-ray  tube  movement.  However, 
with  the  partial  isocentric  motion  of  the  x-ray  tube  in  breast  tomosynthesis,  small  objects 
such  as  microcalcifications  appear  slightly  blurred  (for  instance,  about  1~  1 0  pixels  in  blur 
50  for  a  microcalcification  in  a  human  breast)  in  traditional  SAA  images  in  the  direction 
perpendicular  to  the  direction  of  tube  motion.  Furthermore,  out-of-plane  objects  manifest 
themselves  as  arc-shaped  blurs  due  to  the  isocentric  motion.  Some  digital  breast 
tomosynthesis  algorithms  reported  in  the  literature  utilize  a  traditional  one-dimensional 
SAA  method  that  is  not  wholly  suitable  for  isocentric  motion.  In  this  paper,  a  point-by- 
55  point  back  projection  (BP)  method  is  described  and  compared  with  traditional  SAA  for 
the  important  clinical  task  of  evaluating  morphology  of  small  objects  such  as 
microcalcifications.  Impulse  responses  at  different  3-D  locations  with  five  different 
combinations  of  imaging  acquisition  parameters  were  investigated.  Reconstruction 
images  of  microcalcifications  in  a  human  subject  were  also  evaluated.  Results  showed 
60  that  with  traditional  SAA  and  50°  angular  range  of  tube  movement,  the  in-plane  blur  and 
arc-shaped  out-of-plane  artifacts  were  obvious  for  objects  farther  away  from  x-ray 
source.  In  a  human  subject,  the  appearance  of  calcifications  was  blurred  in  the  direction 
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orthogonal  to  the  tube  motion  and  the  out-of-plane  artifact  of  calcifications  was 
curvilinear  with  traditional  SAA.  With  point-by-point  BP,  the  appearance  of 
65  calcifications  was  sharper.  The  point-by-point  BP  method  demonstrated  improved 
rendition  of  microcalcifications  in  the  direction  perpendicular  to  the  tube  motion 
direction.  With  wide  angles  or  for  imaging  of  larger  breasts,  this  point-by-point  BP  rather 
than  the  traditional  SAA  should  also  be  considered  as  the  basis  of  further  deblurring 
algorithms  that  work  in  conjunction  with  the  BP  method. 

70 

KEYWORDS 

mammography,  tomosynthesis,  3D  reconstruction,  Shift-And-Add  (SAA),  Back 
Projection  (BP),  microcalcifications 
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1.  INTRODUCTION 

75  Breast  cancer  is  the  most  common  cancer  among  women.  Currently,  mammography  is 
the  most  important  and  efficacious  tool  for  the  early  detection  of  breast  cancer.1 
However,  limitations  of  mammography  have  been  well  publicized,  such  as  20%  false 
negative  rate  2’3,  many  callbacks  from  screening,  and  low  positive  predictive  value  of 
about  15-34%  from  biopsy.4’5  It  can  be  difficult  for  conventional  two-dimensional 
80  mammography  to  distinguish  a  cancer  from  overlying  breast  tissues. 

Digital  breast  tomosynthesis  is  a  three-dimensional  imaging  technique  that  provides  an 
arbitrary  set  of  reconstruction  planes  in  the  breast  from  a  limited-angle  series  of 
projection  images  when  the  x-ray  tube  moves  6.  There  are  a  variety  of  tomosynthesis 
85  reconstruction  algorithms,  including  the  image-stretching  method  proposed  by  Niklason 
and  colleagues  ,  Maximum  Likelihood  iterative  algorithm  by  Wu  et  al.  ,  tuned-aperture 
computed  tomography  (TACT)  reconstruction  methods  developed  by  Webber  and 
investigated  by  Suryanarayanan  et  al10,  n,  algebraic  reconstruction  techniques  (ART)  29 ’ 

’  ,  filtered  back  projection  (FBP)  ’  ,  and  matrix  inversion  tomosynthesis  (MITS)  . 

90  Some  of  these  algorithms  depend  on  a  traditional  Shift-And-Add  (SAA)  method  that  is 
appropriate  for  parallel-path  geometries.  For  example,  Niklason  and  colleagues  modified 
the  traditional  shift-and-add  technique  for  mammography  to  stretch  the  image  along  the 
direction  of  x-ray  tube  motion  to  account  for  the  effects  of  magnification  variation  with 
angle,  but  the  correction  necessary  along  the  direction  perpendicular  to  the  tube  motion 
95  was  not  taken  into  account7.  Suryanarayanan  et  al  applied  Webber’s  TACT  method  to 
breast  tomosynthesis  reconstruction  l0,  n,  and  used  traditional  SAA.  The  MITS  technique 
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developed  in  our  laboratory  has  been  investigated  for  breast  tomosynthesis  using  a 

13 

traditional  SAA  algorithm  as  the  basis  for  subsequent  matrix  inversion  deblurring. 

100  Traditional  SAA  is  appropriate  for  parallel-path  tube  movement  when  the  path  of  the 
tube  lies  in  a  plane  parallel  to  the  plane  of  the  detector  6.  However,  the  partial  isocentric 
motion  of  the  tube  in  breast  tomosynthesis  causes  a  non-parallel  motion.  While  the 
effects  due  to  isocentric  motion  are  small  for  most  objects,  the  use  of  SAA  methods 
introduces  morphological  distortions  with  small  objects  such  as  microcalcifications. 

105  Therefore,  a  simple  SAA  reconstruction  algorithm  is  not  entirely  suitable  for  breast 
tomosynthesis.  This  issue  of  non-parallel  motion  is  addressed  in  point-by-point  Back 
Projection  (BP)  methodologies  9,  but  its  impact  has  largely  not  been  evaluated  with 
algorithms  that  rely  on  simply  traditional  SAA  approaches.  This  paper  demonstrates  the 
importance  of  point-by-point  corrections  for  isocentric  motion  in  digital  breast 

110  tomosynthesis  by  examining  how  the  morphology  of  microcalcification  reconstructions 
changes  relative  to  a  traditional  SAA  method  that  does  not  employ  point-by-point 
corrections. 

In  this  paper,  a  point-by-point  BP  correction  method  is  described  and  compared  with 

115  traditional  SAA  by  analysis  of  impulse  response.  Impulse  responses  at  different  3-D 
locations  with  five  different  combinations  of  imaging  acquisition  parameters  were 
investigated.  In  addition,  reconstructed  images  of  a  calcification  in  a  human  subject  were 
evaluated  to  demonstrate  the  improvement  in  the  morphology  of  microcalcifications 
associated  with  the  point-by-point  BP  correction  method. 

120 
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120  2.  METHODS 

A.  Breast  tomosynthesis  system 

A  selenium-based  direct  conversion  Siemens  Mammomat  Novation  DR  prototype  system 
was  modified  to  be  used  as  the  breast  tomosynthesis  acquisition  system."  The  detector 
area  was  24cmx  30cm  (2816  x  3584  pixels),  with  a  pixel  pitch  of  85  pm  (different  from 
125  that  used  in  the  clinical  digital  mammography  system  from  the  same  manufacturer).  The 
exposure  and  readout  cycle  was  0.8  sec  per  image.  Several  different  modes  were 
provided  to  choose  from  different  available  projection  numbers,  total  angular  range  and 
speeds.  Figure  1  shows  a  diagram  of  the  breast  tomosynthesis  imaging  system.  During 
the  tomosynthesis  procedure,  the  x-ray  tube  moves  automatically  along  an  arc  above  the 
130  chest  wall  to  acquire  up  to  49  projection  images  with  a  total  angular  range  of  0-50°  at  the 
rotation  center.  A  continuous  x-ray  motion  was  employed.  The  rotation  center  to  detector 
distance  R  is  6  cm.  A  compression  paddle  is  used  to  keep  the  object  still. 

B.  Traditional  Shift- And- Add  (SAA)  algorithm 

6  12 

135  The  traditional  Shift- And-Add  (SAA)  tomosynthesis  reconstruction  algorithm  ’  "  is  a 
common  mathematical  method  to  line  up  each  projection  image  based  on  its  relative  shift 
to  generate  reconstruction  slices  at  specified  depths.  When  the  x-ray  tube  moves,  objects 
at  different  heights  above  the  detector  will  be  projected  onto  the  detector  at  positions 
depending  on  the  relative  heights  of  the  objects. 

140 

In  order  to  reconstruct  slice  images  of  the  breast,  each  projection  image  should  be  shifted 
by  an  amount  appropriate  for  the  plane  of  reconstruction.  If  the  detector  remains 
stationary  and  the  tube  moves  in  a  plane  that  is  parallel  to  the  detector  plane,  the 
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magnification  of  objects  depends  only  on  the  height  of  the  object.  With  the  traditional 
145  Shift-And-Add  (SAA)  algorithm  for  breast  tomosynthesis  reconstruction,  shift  amounts 
for  each  projection  plane  are  calculated  one  dimensionally  along  the  axis  of  x-ray  tube 
movement.  In  this  paper,  the  shift  amount  was  calculated  based  on  projected  positions 
from  central  points  of  each  reconstruction  plane.  The  shifted  planes  were  added  together 
to  emphasize  structures  in  the  in-focus  plane  and  blur  out  structures  in  other  planes. 

150 

As  shown  in  figure  1,  plane  S  represents  a  reconstruction  plane  at  a  height  of  Z  above  the 
detector  surface.  When  the  x-ray  tube  moves,  objects  in  plane  S  will  be  projected  onto  the 
detector  surface.  For  a  specific  projection  image  from  angle  0,  in  order  to  shift  the 
projection  image  to  line  up  structures  from  plane  S,  the  one-dimensional  SAA  algorithm 
155  uses  the  shift  amount  calculated  as: 

Mft,(Z)  =  P,(Z)  =  LS  inO- - /  (1) 

L  ■  cos  0  +  (R  —  Z) 

where  L  is  length  of  the  rotation  arm,  and  R  is  the  height  of  the  rotation  axis  from  the 
detector  surface.  One  can  obtain  the  reconstruction  plane  S  as  the  average  of  all  N  shifted 
projection  images: 

160  T(x,y)  =  -*-£/,  (x,y)  0  S[y  -  shift,  (Z)]  (2) 

A  i= i 

C.  Point-by-point  Back  Projection  (BP)  algorithm 

Because  of  the  isocentric  motion  of  the  x-ray  tube,  a  shift  actually  occurs  in  both  x  and  y 
165  directions  on  each  projection  image.  Figure  2  shows  the  arc  path  of  motion  when  the  x- 
ray  tube  moves  along  the  y  axis.  Point  A  represents  a  single  structure  on  a  certain 
reconstruction  plane.  P„  Pj,  Pk  are  the  actual  projected  locations  of  point  A  on  the  detector 
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with  different  x-ray  tube  locations  of  Th  Tj,  Tk.  The  actual  path  of  projected  locations 
follows  a  two-dimensional  arc  rather  than  a  one-dimensional  line.  Therefore,  to 
170  reconstruct  a  single  pixel  on  a  reconstruction  plane  at  certain  height  about  the  detector, 
the  shift  amount  should  be  considered  along  both  x  and  y  directions. 


With  the  point-by-point  BP  algorithm,  shift  amounts  for  every  pixel  location  on  each 
reconstructed  plane  were  computed,  taking  into  account  the  2D  arc  projection  location  of 
175  reconstructed  objects  in  each  plane.  In  figure  2,  ( Ax,Ay,Az )  represent  coordinates  of 

point  A.  (Txi,Ty t,Tz.)  represents  the  tube  position  along  the  x,  y,  z  axes  when  the  tube 
moves  to  position  Tj.  ( Pxi,Pyj,Pzi )  represents  projected  coordinates  of  point  A  on  the 
projection  image.  One  can  calculate  the  two-dimensional  shift  amount  as: 


px,  =  Tx,  - 


P)’,  =  Tvt  - 


(TXj  -Ax)(Tzi  -PZj) 
Tz .  -  Az 

C Ty I  -  Ay)(Tz ,.  -  Pzi ) 
Tz ..  -  Az 


(3) 


180  Since  P,  is  located  on  the  detector,  one  can  define  Pz,  =  0 .  Thus,  the  above  formula  can 


be  simplified  as: 


Pxt  =  Txx  + 


Py,  =  Ty,  + 


( Tx j  -  Ax)Tzl 
Tzx  -  Az 

( Ty i  -  Ay)Tzi 

Tz,  -  Az 


(4) 


The  final  pixel  value  of  point  A  in  the  tomosynthesized  reconstruction  was  calculated 
1  N 

as:  —  X7(^), where  7(^)  is  the  pixel  value  at  a  given  location  on  the  /  projection 

N  i=1 

185  image,  and  N  is  the  total  number  of  projection  images.  In  this  paper,  bilinear  interpolation 
was  used  to  address  the  issue  of  partial  pixel  locations.  Computation  times  for  the  point- 
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by-point  BP  algorithm  are  roughly  comparable  to  the  SAA  method.  With  a  computer  of 
800  MHz  CPU  and  UNIX  operating  system,  it  takes  less  than  5  minutes  for  either  the 
point-by-point  BP  or  traditional  SAA  reconstruction. 

190 

D.  Impulse  response  analysis 

A  single  delta  function  was  simulated  by  ray-tracing  method  as  the  input  impulse  to 
investigate  the  sharpness  of  reconstructed  in-plane  structures  and  to  see  how  the 
traditional  SAA  and  point-by-point  BP  algorithm  differ  from  each  other.  Two  different 
195  impulse  locations  were  investigated  in  this  paper:  1)  an  impulse  that  is  exactly  underneath 
the  x-ray  source  (near  the  chest  wall)  and  in  a  defined  reconstructed  plane  (40mm  above 
the  detector  surface);  2)  an  impulse  that  is  approximately  4  cm  away  from  the  chest  wall 
and  in  a  defined  reconstructed  plane  (40mm  above  the  detector  surface).  Parameters  of 
the  digital  breast  tomosynthesis  device  described  in  section  2A  were  used  for  geometries 
200  of  the  simulation.  Five  different  combinations  of  acquisition  parameters  including 
projection  image  numbers  and  total  angular  range  were  applied:  1)  13  projections  with 
25°  angular  range;  2)  13  projections  with  50°  angular  range;  3)  25  projections  with  25° 
angular  range;  4)  25  projections  with  50°  angular  range;  and  5)  49  projections  with  50° 
angular  range.  For  each  impulse  location  and  combination  of  acquisition  parameters,  two 
205  datasets  of  projection  images  were  simulated  by  ray-tracing  method:  1)  background-only: 
only  1/r2  shading  difference  for  each  pixel  on  projection  images  was  taken  into  account  (r 

is  the  distance  from  the  x-ray  source  to  each  pixel  location);  2)  impulse-added:  projection 

2 

images  with  simulated  impulse  and  the  1/r  shading  difference.  Other  system  blur  and 
noise  issues  were  not  addressed  in  this  paper  to  focus  on  the  contribution  of  the  blur  due 
210  to  the  isocentric  motion.  During  ray-tracing,  if  the  simulated  impulse  was  projected  onto 
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non-integer  pixel  location  on  the  detector  surface,  linear  interpolation  of  the  projected 
impulse  among  four  neighboring  pixels  was  performed. 

Traditional  SAA  and  point-by-point  BP  reconstruction  algorithms  were  applied  to  both 
215  impulse-added  and  background-only  simulated  tomosynthesis  projection  sequences.  A 
reconstruction  plane  spacing  of  1  mm  was  used.  Reconstruction  images  from 
background-only  projections  were  subtracted  from  reconstructions  of  impulse-added 
projections  to  eliminate  background  shading  effects.  The  impulse  responses  were 
nonnalized  based  on  the  ideal  condition  when  the  impulse  is  exactly  located  underneath 
220  the  x-ray  source. 

E.  Human  subject  images 

Human  subject  images  have  been  acquired  on  our  prototype  breast  tomosynthesis  system 
under  an  IRB-approved  protocol.  Images  of  one  human  subject  with  notable 
225  calcifications  were  reconstructed  with  the  traditional  SAA  and  point-by-point  BP 
methods  to  demonstrate  the  effect  of  the  point-by-point  BP  method  on  reconstructed 
calcification  morphological  appearance.  A  tomosynthesis  sequence  was  acquired  with 
twenty-five  projection  images  and  a  total  angular  range  of  50°.  The  radiographic 
technique  for  breast  tomosynthesis  was  selected  using  technique  optimization  procedures 
230  reported  previously  .  The  target/filter  for  tomosynthesis  exams  is  always  W/Rh.  The 
kVp  is  selected  to  maximize  a  figure  of  merit  (signal  difference  to  noise  ratio  squared  per 
unit  dose)  for  a  given  breast  thickness  and  density.  Then  the  mAs  can  be  chosen  to 
maintain  image  quality  while  reducing  dose  when  compared  against  the  conventional 
Mo/Mo  or  Mo/Rh  technique,  or  alternatively  to  improve  quality  while  maintaining  dose. 
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235  For  this  specific  subject,  we  chose  to  do  the  latter.  By  using  the  tomosynthesis  technique 
of  W/Rh  at  28  kVp  (HVL  0.50  mm  Al)  and  1 12  mAs  for  this  100%  fatty,  33  mm  breast, 
we  maintained  the  same  dose  as  the  conventional  left  cranio  caudal  (LCC)  mammogram. 
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240  3.  RESULTS 

A.  Impulse  responses 

Figures  3,  4,  and  5  show  the  impulse  response  results  with  a  simulated  impulse  located  in 
a  defined  reconstruction  plane  at  40  mm  above  the  detector,  and  approximately  4  cm 
away  from  the  chest  wall.  Figure  3  shows  results  from  simulated  acquisition  parameters 
245  of  25  projections  and  50°  total  tube  angular  movement.  Figure  4  shows  results  from  13 
projections  and  50°  angular  movement.  Figures  5  shows  results  from  25  projections  and  a 
narrower  angular  range  of  25°.  On  figures  3  through  5,  (a)  and  (c)  give  corresponding 
values  in  a  plane  at  the  exact  height  of  the  impulse’s  location;  (b)  and  (d)  give  the 
impulse  response  of  reconstruction  planes  1  mm  lower  than  the  impulse’s  location;  (a) 
250  and  (b)  are  results  from  point-by-point  BP,  (c)  and  (d)  are  results  from  traditional  SAA. 
The  x  and  v  axes  give  the  pixel  location  on  the  reconstruction  plane,  and  the  plot  displays 
the  normalized  amplitude  of  the  response.  X  axis  represents  the  direction  orthogonal  to 
tube’s  motion  direction.  Y  axis  represents  the  tube’s  motion  direction.  Only  a  40x40 
pixel  region  close  to  the  impulse  is  shown  for  clarity. 

255 

One  can  see  that  with  traditional  SAA,  when  the  impulse  is  4  cm  away  from  the  chest 

wall  and  the  x-ray  source  moves  along  a  wider  angular  range  of  50°,  the  in-plane 

response  is  noticeably  blurred  and  multiple  peaks  exist  in  a  direction  that  is  perpendicular 

to  the  direction  of  tube  movement  (Fig.  3c  and  4c),  reflecting  the  uncorrected  partial 

260  isocentric  tube  motion.  With  point-by-point  BP,  the  in-plane  response  is  much  sharper 

(Fig.  3a  and  4a),  and  the  out-of-plane  blur  responses  are  less  curved  compared  with  that 

of  traditional  SAA  (Fig.  3b,  2d  and  4b,  4d).  When  the  number  of  projection  images 

decreases  to  13,  the  in-plane  response  and  out-of-plane  blur  become  discrete  (Fig.  4b,  4d) 
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due  to  limited  projection  numbers.  With  a  narrower  angular  range  of  25°,  the  differences 
265  between  traditional  SAA  and  point-by-point  BP  are  less  obvious  (Fig.  5).  However,  one 
can  still  say  that  the  in-plane  response  of  point-by-point  BP  is  higher  and  sharper  than 
that  of  traditional  SAA. 

Table  1  gives  the  full  width  at  a  half  maximum  (FWHM)  measurement  and  full  width  at  a 
270  tenth  maximum  values  (FWTM)  of  the  in-plane  impulse  responses  along  two  orthogonal 
directions  when  the  impulse  is  located  at  40  mm  above  the  detector  surface  and  near 
chest  wall.  Table  2  gives  the  same  measurements  of  FWHM  and  FWTM  when  the 
impulse  is  located  4  cm  away  chest  wall.  When  the  impulse  is  located  near  the  chest  wall 
(underneath  the  x-ray  tube),  there  is  only  a  small  difference  (<  1  pixel  size)  in  FWHM 
275  and  FWTM  values  between  traditional  SAA  and  point-by-point  BP  for  each  combination 
of  acquisition  parameters.  With  a  narrower  angle  of  25°,  differences  are  small  too  (less 
than  1  pixel  size).  However,  with  a  wider  tube  angular  movement  range  of  50°,  when  the 
impulse  is  located  4  cm  away  the  chest  wall,  major  differences  exist  along  the  X  axis 
(direction  orthogonal  to  tube’s  motion  direction).  Multiple  peaks  and  blurs  appear  along 
280  this  direction  (Fig.  3c  and  4c).  Due  to  multiple  peaks  along  the  X  axis,  FWHM  cannot 
adequately  represent  the  real  extension  of  impulse  response;  for  that  reason,  the 
measurement  of  FWTM  was  also  provided.  The  difference  in  FWTM  values  is  as  large  as 
9  pixels  between  traditional  SAA  and  point-by-point  BP. 

285  B.  Human  subject  images 

Figures  6  shows  regions  of  interest  containing  two  calcifications  from  images  of  the  left 
(L)  breast  of  the  human  subject.  Reconstructed  structures  by  traditional  SAA  and  point- 
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by-point  BP  methods  are  compared.  A  12. 75x  12.75  mm  region-of-interest  (ROI)  is 
demonstrated.  Compared  with  traditional  SAA,  edges  of  calcifications  are  sharper  and 
290  better  focused  with  point-by-point  BP  (Fig.  6b  and  6e).  With  traditional  SAA,  out-of¬ 
plane  structures  appeared  curvilinear  and  indistinct  (Fig.  6a  and  6c).  With  point-by-point 
BP,  the  in-plane  structures  are  sharper,  and  the  out-of-plane  structures  appear  linear  and 
curvilinear  (Fig.  6d  and  6f).  A  quantitative  measurement  of  the  shape  and  width  of  the 
punctate  calcification  (left-most  calcification)  on  this  ROI  is  given  in  table  3.  Y  axis 
295  corresponds  to  the  direction  of  x-ray  tube  motion,  and  X  axis  is  the  direction  orthogonal 
to  the  tube  motion  direction.  One  can  see  that,  compared  with  the  traditional  SAA,  point- 
by-point  BP  provided  clearer  edge  shape  and  narrower  width  along  X  axis  that  is 
perpendicular  to  the  x-ray  tube  motion  direction. 

300  Figure  7a  is  a  low  dose  middle  projection  image  of  the  same  human  breast  when  the  x-ray 
tube  was  positioned  at  the  0°  position.  Figures  7b  and  7c  are  reconstructed  slice  images 
by  traditional  one-dimensional  SAA  and  point-by-point  BP  respectively,  at  a  height  of 
7.5  mm  above  the  detector.  One  can  find  that  the  rendition  of  a  solitary  calcification 
(when  out  of  the  reconstructed  plane)  appears  more  curvilinear  in  the  SAA  image 
305  compared  with  that  in  the  point-by-point  BP  image.  Also,  the  reconstructed  size  of  the 
breast  by  traditional  one-dimensional  SAA  is  larger  than  that  from  point-by-point  BP 
(Fig.  7b,  7c).  This  is  due  to  the  uncorrected  magnification  difference  with  traditional 
SAA  for  structures  at  different  locations  in  the  reconstruction  plane.  When  the  structures 
are  located  at  difference  heights  above  the  detector  surface,  differences  in  magnification 
310  exist.  Structures  at  higher  locations  above  the  detector  surface  will  be  projected  onto  the 
detector  with  a  larger  size.  Traditional  SAA  doesn’t  take  this  magnification  difference 
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into  account.  After  shift-and-add,  structures  at  higher  locations  appear  larger  compared 
with  the  same  size  structure  at  lower  location.  The  reconstructed  structure  size  from 
traditional  SAA  reconstruction  cannot  reflect  the  real  structure  size.  Point-by-point  BP 
315  correctly  addressed  this  magnification  difference  issue  by  calculating  shift  amounts  for 
every  pixel  location  on  each  reconstruction  plane.  Therefore,  the  point-by-point  BP 
reconstructed  plane  in  figure  7c  reflects  the  real  size  and  is  smaller  than  that  from 
traditional  SAA. 

320 


26 


16 


320  4.  DISCUSSION 

With  the  traditional  SAA  algorithm  for  digital  breast  tomosynthesis  reconstruction,  shift 
amounts  for  each  projection  plane  are  calculated  only  along  the  axis  of  the  x-ray  tube’s 
movement.  However,  due  to  the  isocentric  motion  of  x-ray  tube,  the  track  of  projected 
impulse  locations  actually  occurs  in  two  dimensions  on  the  detector.  As  a  result,  the  in- 
325  plane  structures  are  blurred  and  the  out-of-plane  structures  have  curve-shaped 
appearance.  Illustrations  with  impulse  responses  and  reconstructed  human  subject  images 
demonstrated  that  this  is  an  inherent  problem  of  traditional  one-dimensional  SAA  with 
breast  tomosynthesis.  This  problem  is  more  obvious  when  the  object  is  farther  away  from 
the  chest  wall  and  higher  above  the  detector.  Thus,  with  traditional  SAA,  objects  such  as 
330  microcalcifications  appear  slightly  blurred  in  the  direction  perpendicular  to  the  direction 
of  tube  motion,  and  their  apparent  morphology  changes. 

With  point-by-point  BP,  the  artifacts  coming  from  the  isocentric  x-ray  tube’s  movement 
are  corrected.  The  in-plane  structures  are  sharper.  While  reconstructions  of  gross 
335  anatomy  were  adequate  with  either  algorithm,  the  morphology  of  small  structures  such  as 
microcalcifications  reconstruction  requires  a  point-by-point  correction.  Morphological 
artifacts  were  reduced  with  the  point-by-point  correction,  and  rendition  of  small  objects 
such  as  microcalcifications  were  greatly  improved.  These  results  demonstrate  the 
importance  of  using  a  point-by-point  correction  to  remove  isocentric  motion  artifacts  in 
340  tomosynthesis  imaging  of  the  breast,  where  the  morphology  of  microcalcifications  has  an 
important  bearing  on  clinical  decision  making. 
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The  source  of  the  difficulties  with  traditional  SAA  is  the  variable  magnification 
introduced  by  the  partial  isocentric  motion  of  the  x-ray  tube.  The  magnification  of  objects 
345  in  different  reconstruction  slices  varies  from  plane-to-plane  as  a  function  of  the  height  of 
slices  above  the  detector.  Even  within  the  same  reconstruction  plane,  the  magnification 
also  changes  at  different  pixel  locations  due  to  the  partial  isocentric  motion.  Traditional 
SAA  does  not  take  this  issue  into  account.  With  point-by-point  BP,  the  shift  amount  is 
calculated  according  to  the  exact  location  of  the  pixel  in  the  reconstruction  slices,  thereby 
350  addressing  this  issue  of  variable  magnification. 

None  of  methods  commonly  used  for  breast  tomosynthesis  is  truly  spatially  invariant  due 
to  the  partial  isocentric  tube  motion,  although  for  structures  close  to  the  chest  wall  and 
near  the  detector,  the  imaging  system  is  approximately  spatially  invariant.  Therefore, 
355  linear  deblurring  techniques  will  demonstrate  less  uniform  behavior  in  breast 
tomosynthesis  geometries  than  in  parallel-path  tomosynthesis  geometries.  We  found  that 
when  the  impulse  is  located  near  the  chest  wall  or  when  the  total  angular  range  is 
narrower  such  as  25°,  there  is  no  big  difference  (<1  pixel  size)  in  FWHM  and  FWTM 
values  between  traditional  SAA  and  point-by-point  BP  for  a  reconstruction  with  a  height 
360  of  40  mm  above  the  detector.  However,  with  a  wider  angle  of  50°,  even  moving  the 
impulse  4  cm  away  chest  wall  shows  a  difference  (about  9  pixels)  along  the  direction 
orthogonal  to  the  direction  of  tube  motion.  Therefore,  with  a  narrow  angle,  or  for  small  or 
thin  breasts,  the  SAA  algorithm  may  be  tolerated.  However,  with  a  wide  angle  or  large 
breast  size,  the  point-by-point  BP  algorithm  rather  than  traditional  SAA  should  be  used  to 
365  minimize  issues  related  to  isocentric  motion. 
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Deblurring  algorithms  such  as  FBP  or  MITS,  which  are  an  important  component  of  high- 
quality  tomosynthesis  reconstruction,  should  use  the  point-by-point  BP  method  rather 
than  the  SAA  method  to  generate  the  constituent  images  prior  to  deblurring  under  the 
370  same  conditions  of  wide  tube  angle  or  large  breast  size. 
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5.  CONCLUSIONS 

This  work  demonstrates  that  point-by-point  BP  is  an  effective  method  to  reconstruct  3D 
375  tomosynthesis  images  of  the  breast  with  improved  rendition  of  small  structures  such  as 
microcalcifications.  Compared  with  the  traditional  SAA  algorithm,  the  method  of  point- 
by-point  BP  takes  into  account  the  variable  magnification  and  shift  occurring  along  the 
direction  orthogonal  to  tube  movement  due  to  the  isocentric  tube  motion.  Point-by-point 
BP  improves  the  sharpness  and  morphology  of  structures  especially  for  small  objects 
380  such  as  calcifications.  This  may  prove  helpful  to  radiologists  in  discriminating  malignant 
from  benign  microcalcification  patterns,  and  thereby  improve  the  accuracy  of  breast 
cancer  detection. 
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390  Figure  captions 

FIGURE  1.  Breast  Tomosynthesis  Imaging  System  ( O  is  the  rotation  center.  R  is  the 
rotation  center  to  detector  distance.  L  is  the  rotation  arm  length.  Z  is  the  height  of  plane  S 
above  the  detector.) 

395  FIGURE  2.  Two-dimensional  arc  path  from  isocentric  tube  motion 

FIGURE  3.  Traditional  SAA  and  point-by-point  BP  impulse  responses:  25  projection 
images  and  50°  angular  range,  with  the  impulse  located  40mm  above  the  detector  and 
about  4cm  away  from  the  chest  wall,  (a)  and  (b)  give  the  impulse  response  of  point-by- 
400  point  BP;  (c)  and  (d)  give  the  impulse  response  of  traditional  SAA.  (a)  and  (c)  are 
impulse  responses  in  the  plane  at  the  exact  height  of  the  impulse’s  location;  (b)  and  (d) 
give  corresponding  values  in  a  reconstruction  plane  1  mm  below  the  impulse’s  location. 

FIGURE  4.  Traditional  SAA  and  point-by-point  BP  impulse  responses:  13  projection 
405  images  and  50°  angular  range,  with  the  impulse  located  40mm  above  the  detector  and 
about  4cm  away  from  the  chest  wall,  (a)  and  (b)  give  the  impulse  response  of  point-by¬ 
point  BP;  (c)  and  (d)  give  the  impulse  response  of  traditional  SAA.  (a)  and  (c)  are 
impulse  responses  in  the  plane  at  the  exact  height  of  the  impulse’s  location;  (b)  and  (d) 
give  corresponding  values  in  a  reconstruction  plane  1  mm  below  the  impulse’s  location. 
410 

FIGURE  5.  Traditional  SAA  and  point-by-point  BP  impulse  responses:  25  projection 
images  and  25°  angular  range,  with  the  impulse  located  40mm  above  the  detector  and 
about  4cm  away  from  the  chest  wall,  (a)  and  (b)  give  the  impulse  response  of  point-by- 
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point  BP;  (c)  and  (d)  give  the  impulse  response  of  traditional  SAA.  (a)  and  (c)  are 
415  impulse  responses  in  the  plane  at  the  exact  height  of  the  impulse’s  location;  (b)  and  (d) 
give  corresponding  values  in  a  reconstruction  plane  1  mm  below  the  impulse’s  location. 

TABLE  1 .  The  impulse  is  located  near  chest  wall.  Full  width  at  half  maximum  (FWHM) 
and  full  width  at  tenth  maximum  (FWTM)  measurements  of  in-plane  impulse  response 
420  along  two  directions:  tube’s  motion  direction  (Y  axis)  and  direction  orthogonal  to  tube’s 
motion  (X  axis). 

Table  2.  The  impulse  is  located  4  cm  away  chest  wall.  Full  width  at  half  maximum 
(FWHM)  and  full  width  at  tenth  maximum  (FWTM)  measurements  of  in-plane  impulse 
425  response  along  two  directions:  tube’s  motion  direction  (Y  axis)  and  direction  orthogonal 
to  tube’s  motion  (X  axis). 

FIGURE  6.  Reconstructed  12. 75x  12.75mm  ROI  of  a  human  breast  containing 
calcifications,  Z=18  mm  represents  the  plane  height  closet  to  the  location  of  the 
430  calcification:  (a)  Traditional  SAA,  Z=16.5mm;  (b)  Traditional  SAA,  Z=18mm;  (c) 
Traditional  SAA,  Z=19.5mm;  (d)  Point-by-point  BP,  Z=16.5mm;  (e)  Point-by-point  BP, 
Z=1 8mm;  (f)  Point-by-point  BP,  Z=1 9.5mm 

TABLE  3.  Quantitative  measurement  of  the  leftmost  calcification  depicted  in  Fig.  6. 

435 

FIGURE  7.  A  human  breast  demonstrating  a  solitary  calcification:  (a)  Low  dose  middle 

(0°)  projection  image  of  the  tomosynthesis  sequence.  The  spectrum  used  for  the 

tomosynthesis  sequence  was  28  kVp  with  W/Rh  targer/filter  and  1 12.5  mAs  for  a  total  of 
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25  projection  images  and  50°  angular  range,  (b)  Traditional  SAA  reconstructed  slice 
image:  Z=7.5mm.  (c)  Point-by-point  BP  reconstructed  slice  image:  Z=7.5mm 
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FIG.  1.  Breast  Tomosynthesis  Imaging  System  ( O  is  the  rotation  center.  R  is  the  rotation 
center  to  detector  distance.  L  is  the  rotation  ann  length.  Z  is  the  height  of  plane  S  above 
455  the  detector.) 
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FIG.  2.  Two-dimensional  arc  path  from  isocentric  tube  motion 
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FIG.  3.  Traditional  SAA  and  point-by-point  BP  impulse  responses:  25  projection  images 
and  50°  angular  range,  with  the  impulse  located  40mm  above  the  detector  and  about  4cm 
away  from  the  chest  wall,  (a)  and  (b)  give  the  impulse  response  of  point-by-point  BP;  (c) 
465  and  (d)  give  the  impulse  response  of  traditional  SAA.  (a)  and  (c)  are  impulse  responses  in 
the  plane  at  the  exact  height  of  the  impulse’s  location;  (b)  and  (d)  give  corresponding 
values  in  a  reconstruction  plane  1mm  below  the  impulse’s  location. 
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470 

FIG.  4.  Traditional  SAA  and  point-by-point  BP  impulse  responses:  13  projection  images 
and  50°  angular  range,  with  the  impulse  located  40mm  above  the  detector  and  about  4cm 
away  from  the  chest  wall,  (a)  and  (b)  give  the  impulse  response  of  point-by-point  BP;  (c) 
and  (d)  give  the  impulse  response  of  traditional  SAA.  (a)  and  (c)  are  impulse  responses  in 
475  the  plane  at  the  exact  height  of  the  impulse’s  location;  (b)  and  (d)  give  corresponding 
values  in  a  reconstruction  plane  1  mm  below  the  impulse’s  location. 
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FIG.  5.  Traditional  SAA  and  point-by-point  BP  impulse  responses:  25  projection  images 
480  and  25°  angular  range,  with  the  impulse  located  40mm  above  the  detector  and  about  4cm 
away  from  the  chest  wall,  (a)  and  (b)  give  the  impulse  response  of  point-by-point  BP;  (c) 
and  (d)  give  the  impulse  response  of  traditional  SAA.  (a)  and  (c)  are  impulse  responses  in 
the  plane  at  the  exact  height  of  the  impulse’s  location;  (b)  and  (d)  give  corresponding 
values  in  a  reconstruction  plane  1  mm  below  the  impulse’s  location. 
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485 


Acquisition 

Parameters 

Full  Width  at  Tenth  Maximum 
(in  pixels) 

Full  Width  at  Ffalf  Maximum 
(in  pixels) 

Projections 

Total 

angular 

range 

Traditional 

SAA 

Point-by-point 

BP 

Traditional  SAA 

Point-by-point 

BP 

X  axis 

Y  axis 

X  axis 

Y  axis 

X  axis 

Y 

axis 

X  axis 

Y  axis 

13 

25° 

2.8 

2.8 

3.4 

2.3 

1.9 

1.2 

1.7 

1.1 

25 

25° 

2.8 

2.8 

3.4 

2.4 

1.9 

1.2 

1.6 

1.2 

13 

50u 

2.7 

2.7 

3.1 

2.1 

1.4 

1.2 

1.3 

1.2 

25 

50° 

2.7 

2.8 

3.2 

2.3 

1.4 

1.2 

1.4 

1.2 

49 

50° 

2.7 

2.6 

3.2 

2.1 

1.5 

1.2 

1.4 

1.2 

Table  1.  The  impulse  is  located  near  chest  wall.  Full  width  at  half  maximum  (FWHM) 
and  full  width  at  tenth  maximum  (FWTM)  measurements  of  in-plane  impulse  response 
490  along  two  directions:  tube’s  motion  direction  (Y  axis)  and  direction  orthogonal  to  tube’s 
motion  (X  axis). 
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Acquisition 

Parameters 

Full  Width  at  Tenth  Maximum 
(in  pixels) 

Full  Width  at  Half  Maximum 
(in  pixels) 

Projections 

Total 

angular 

range 

Traditional 

SAA 

Point-by-point 

BP 

Traditional  SAA 

Point-by-point 

BP 

X  axis 

Y  axis 

X  axis 

Y  axis 

X  axis 

Y 

axis 

X  axis 

Y  axis 

13 

25u 

3.7 

2.8 

2.8 

2.5 

2.0 

1.2 

1.2 
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50° 

11.4 

2.8 

2.6 

2.1 

(~13)* 

Multiple 

peaks 

1.2 

1.2 

1.2 

*  For  the  50°  scan,  multiple  discrete  peaks  in  the  impulse  response  make  FWHM  not 


meaningful  as  a  descriptor  of  impulse  width.  For  these  cases,  the  approximate  overall 
495  width  of  the  impulse  is  provided. 


Table  2.  The  impulse  is  located  4  cm  away  chest  wall.  Full  width  at  half  maximum 
(FWHM)  and  full  width  at  tenth  maximum  (FWTM)  measurements  of  in-plane  impulse 
500  response  along  two  directions:  tube’s  motion  direction  (Y  axis)  and  direction  orthogonal 
to  tube’s  motion  (X  axis). 
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FIG. 6.  Reconstructed  12. 75x  12.75mm  ROI  of  a  human  breast  containing  calcifications, 
Z=18  mm  represents  the  plane  height  closet  to  the  location  of  the  calcification:  (a) 
Traditional  SAA,  Z=  16.5mm;  (b)  Traditional  SAA,  Z=18mm;  (c)  Traditional  SAA, 
510  Z=19.5mm;  (d)  Point-by-point  BP,  Z=16.5mm;  (e)  Point-by-point  BP,  Z=18mm;  (f) 

Point-by-point  BP,  Z=  19.5mm 
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Z=  16.5mm 

Z=1 8  mm 

Z=1 9.5mm 

Shape 

X 

width 

(in 

pixels) 

Y 

width 

(in 

pixels) 

Shape 

X 

width 

(in 

pixels) 

Y 

width 

(in 

pixels) 

Shape 

X 

width 

(in 

pixels) 

Y 

width 

(in 

pixels) 

Traditional 

SAA 

curvilinear 

10 

21 

punctate 

blurred 

9 

11 

indistinct 

11 

16 

Point-by¬ 
point  BP 

linear 

7 

20 

punctate 

sharp 

7 

9 

curvilinear 

7 

17 

Table  3.  Quantitative  measurement  of  the  leftmost  calcification  depicted  in  Fig.  6. 
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(b)  (c) 


FIG.  7.  A  human  breast  demonstrating  a  solitary  calcification:  (a)  Low  dose  middle  (0°) 
projection  image  of  the  tomosynthesis  sequence.  The  spectrum  used  for  the 
tomosynthesis  sequence  was  28  kVp  with  W/Rh  targer/filter  and  1 12.5  mAs  for  a  total  of 
525  25  projection  images  and  50°  angular  range,  (b)  Traditional  SAA  reconstructed  slice 

image:  Z=7.5mm.  (c)  Point-by-point  BP  reconstructed  slice  image:  Z=7.5mm 
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TWO-DIMENSIONAL  SHIFT-AND-ADD  (SAA)  ALGORITHM  FOR  DIGITAL  BREAST 
TOMOSYNTHESIS  RECONSTRUCTION 


Introduction 

Breast  cancer  is  the  most  common  cancer  among  women.  Currently,  mammography  is  the  most  important  and 
efficacious  tool  for  the  early  detection  of  breast  cancer  [7].  However,  limitations  of  mammography  have  been  well 
publicized,  such  as  20%  false  negative  rate  [8'  9\  many  callbacks  from  screening,  and  low  positive  predictive  value  of 
about  15-34%  from  biopsy  [10'  "l  It  is  very  difficult  for  conventional  two-dimensional  mammography  to  distinguish 
a  cancer  from  overlying  breast  tissues. 

Digital  breast  tomosynthesis  is  a  three-dimensional  imaging  technique  that  provides  an  arbitrary  set  of 
reconstruction  planes  in  the  breast  from  limited-angle  series  of  projection  images  when  x-ray  tube  moves  [1"5].  Shift- 
And-Add  (SAA)  tomosynthesis  reconstruction  algorithm  is  a  common  mathematical  method  to  line  up  each 
projection  image  based  on  its  shifting  amount  to  generate  the  reconstruction  slices.  Various  efforts  have  been  made 
on  SAA  reconstruction  for  breast  imaging  [1,2’3’5].  However,  among  those  efforts,  shift  amounts  for  each  projection 
plane  was  calculated  only  along  the  axis  of  x-ray  tube’s  movement.  This  brings  inaccuracy  to  the  reconstructed  slice 
images.  This  project  aims  to  investigate  a  two-dimensional  correction  for  the  SAA  algorithm  to  correctly  calculate 
the  shift  amount  based  on  isocentric  tube  motion,  in  order  to  generate  more  accurate  and  reliable  reconstruction 
images. 

Methods  and  Materials 

A.  Breast  tomosynthesis  system 

A  selenium-based  direct  conversion  Siemens  Mammomat  Novation  prototype  system  with  the  pixel  size  of  85  pm 
was  used  as  the  breast  tomosynthesis  acquisition  system.  It  can  generate  tomosynthesis  sequences  of  up  to  49 
projection  images  with  a  total  angular  range  of  0-50  degrees.  A  few  tomosynthesis  sequences  of  human  subjects 
were  acquired  for  investigation  and  comparison. 

B.  Two-dimensional  SAA  algorithm 

In  order  to  reconstruct  slice  images  of  the  breast,  each  projection  image  should  be  lined  up  based  on  the  shift  amount 
appropriate  for  the  plane  of  reconstruction.  The  shift  amount  of  each  projection  image  was  calculated  according  to 
the  relative  positions  of  the  projection  image  and  geometric  parameters  of  the  tomosynthesis  acquisition  system.  In¬ 
focus  structures  in  a  given  plane  were  lined  up  for  reconstruction. 

During  tomosynthesis  sequence  acquisition,  the  x-ray  tube  moves  above  the  chest  wall  along  an  arc.  Figure  1  shows 
the  arc  path  of  motion  when  x-ray  tube  moves  along  x  axis.  Point  P  represents  a  single  structure  on  a  certain 
reconstruction  plane.  Ph  P2,  P  ;  are  the  actual  projected  locations  of  point  P  on  the  detector  for  different  x-ray  tube 
locations  of  Ab  A 2,  A  3  specifically.  The  actual  path  of  the  projected  location  follows  a  two-dimensional  arc  rather 
than  a  one-dimensional  line.  However,  with  the  traditional  Shift- And-Add  (SAA)  algorithm  for  breast  tomosynthesis 
reconstruction,  shift  amounts  for  each  projection  plane  are  calculated  only  along  one  axis  of  x-ray  tube’s  movement. 
As  a  result,  small  objects  such  as  microcalcifications  appear  slightly  blurred  in  the  direction  perpendicular  to  the 
direction  of  tube  motion. 

Because  of  the  isocentric  motion  of  the  x-ray  tube,  the  shift  amount  occurs  at  both  x  and  y  two-dimensional 
directions  on  each  projection  image.  For  our  two-dimensional  SAA  method,  shift  amounts  for  every  pixel  location 
on  each  reconstruction  plane  were  computed,  taking  into  account  the  2D  arc  projection  location  of  reconstructed 
objects  in  each  plane.  Bilinear  interpolation  was  used  for  partial  pixel  locations. 

C.  Impulse  response  analysis 

Impulses  at  different  3-D  locations  were  simulated  to  investigate  the  sharpness  of  reconstructed  in-plane  structures 
and  the  effectiveness  at  removing  out-of-plane  blur.  Datasets  with  49,  25  and  13  projection  images  of  the  impulse 
were  generated  with  a  total  angular  movement  of  25  and  50  degrees  of  the  simulated  x-ray  point  source.  Four  groups 
of  impulses  were  simulated:  1)  impulse  that  is  exactly  underneath  the  x-ray  source  (near  the  chest  wall)  and  in  a 
defined  reconstructed  plane;  2)  impulse  that  is  exactly  underneath  the  x-ray  source  (near  the  chest  wall)  and  halfway 
between  reconstructed  planes;  3)  impulse  that  is  approximately  4  cm  away  from  the  chest  wall  and  in  a  defined 
reconstructed  plane;  4)  impulse  that  is  approximately  4  cm  away  from  the  chest  wall  and  hallway  between 
reconstructed  planes.  The  impulse  response  in  the  spatial  domain  was  analyzed  for  evaluation. 
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Reconstruction  Plane 
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Figure  1.  Arc  path  of  motion  when  x-ray  tube  moves  along  x  axis 


Results 

According  to  impulse  response  analysis,  two-dimensional  SAA  demonstrated  the  improvement  in  the  direction  that 
is  perpendicular  to  the  tube  motion  direction.  According  to  human  subjects  results,  the  appearance  of  calcifications 
from  two-dimensional  SAA  was  sharper  than  traditional  SAA  at  the  direction  orthogonal  to  the  tube  motion 
direction.  Out-of-plane  artifacts  of  calcifications  changed  from  curved  to  be  straight. 

Figure  2  shows  two  reconstruction  slices  for  a  human  subject.  Figure  3  shows  reconstructed  ROI  images  of  the  same 
human  subject.  The  tomosynthesis  sequence  was  acquired  with  25  projections  and  +/-25°angular  range,  with  a  total 
exposure  of  1.5  times  of  a  single  CC  view.  Figure  3(a),  (b)  are  reconstructed  ROIs  of  a  calcification  on  a  single  slice 
from  traditional  SAA  and  2D  SAA  specifically.  We  can  find  that  the  appearance  of  calcifications  from  two- 
dimensional  SAA  was  sharper  and  more  accurately  defined.  Figure  3(c),  (d)  show  the  out-of-plane  blur  of  a 
calcification.  Out-of-plane  artifacts  of  calcifications  changed  from  curved  to  be  straight  after  2D  SAA. 

We  also  found  that,  when  x-ray  tube  moves  above  the  breast  along  an  arc,  for  different  locations  of  structures  on  a 
single  reconstruction  plane,  the  magnification  difference  existed.  This  is  more  substantial  for  structures  further  away 
from  the  chest  wall.  Two-dimensional  SAA  correctly  addressed  this  issue  by  calculating  shift  amounts  for  every 
pixel  location  on  each  reconstruction  plane. 
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Figure  2(a). 

2D  SAA: 

Single  reconstruction  slice 
with  mass  on  it 


Figure  2(b). 

2D  SAA: 

Single  reconstruction  slice 
with  big  calcification  on  it 
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(a)  (b) 


\ 


(c)  (d) 


Figure  3.  Reconstructed  ROI  of  a  human  subject: 

(a)  Traditional  SAA:  calcification 

(b)  Two-dimensional  SAA:  calcification 

(c)  Traditional  SAA:  out-of-plane  blur  of  calcification 

(d)  Two-dimensional  SAA:  out-of-plane  blur  of  calcification 


Conclusions 

This  work  demonstrates  that  the  two-dimensional  SAA  is  an  effective  method  to  reconstruct  3D  tomosynthesis 
images  of  the  breast.  Compared  with  the  traditional  SAA  algorithm,  the  new  method  corrects  for  two-dimensional 
shift  amounts  coming  from  the  isocentric  x-ray  tube’s  movement.  This  provides  more  accurate  and  reliable  results 
than  other  SAA  algorithms. 

The  two-dimensional  SAA  does  appear  to  improve  the  sharpness  and  morphology  of  calcifications.  However,  it  is 
not  spatially  invariant  method,  therefore  may  be  more  difficult  to  be  used  with  linear  deblurring  techniques. 
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ABSTRACT 

As  a  new  three-dimensional  imaging  technique,  digital  breast  tomosynthesis  allows  the  reconstruction  of  an  arbitrary 
set  of  planes  in  the  breast  from  a  limited-angle  series  of  projection  images.  Though  several  tomosynthesis  algorithms 
have  been  proposed,  no  complete  optimization  and  comparison  of  different  tomosynthesis  acquisition  techniques  for 
available  methods  has  been  conducted  as  of  yet.  This  paper  represents  a  methodology  of  noise-equivalent  quanta 
NEQ  (f)  analysis  to  optimize  and  compare  the  efficacy  of  tomosynthesis  algorithms  and  imaging  acquisition 
techniques  for  digital  breast  tomosynthesis.  It  combines  the  modulation  transfer  function  (MTF)  of  system  signal 
performance  and  the  noise  power  spectrum  (NPS)  of  noise  characteristics.  It  enables  one  to  evaluate  the 
performance  of  different  acquisition  parameters  and  algorithms  for  comparison  and  optimization  purposes.  An 
example  of  this  methodology  was  evaluated  on  a  selenium-based  direct-conversion  flat-panel  Siemens  Mammomat 
Novation  prototype  system.  An  edge  method  was  used  to  measure  the  presampled  MTF  of  the  detector.  The  MTF 
associated  with  the  reconstruction  algorithm  and  specific  acquisition  technique  was  investigated  by  calculating  the 
Fourier  Transform  of  simulated  impulse  responses.  Flat  field  tomosynthesis  projection  sequences  were  acquired  and 
then  reconstructed.  A  mean-subtracted  NPS  on  the  reconstructed  plane  was  studied  to  remove  fixed  pattern  noise. 
An  example  of  the  application  of  this  methodology  was  illustrated  in  this  paper  using  a  point-by-point  Back 
Projection  correction  (BP)  reconstruction  algorithm  and  an  acquisition  technique  of  25  projections  with  25  degrees 
total  angular  tube  movement. 

Keywords:  mammography,  tomosynthesis,  breast  imaging,  modulation  transfer  function  (MTF),  noise  power 
spectrum  (NPS),  Noise-equivalent  quanta  (NEQ),  impulse  response,  Back  Projection  (BP) 


1.  INTRODUCTION 


Digital  breast  tomosynthesis  is  a  new  technique  to  reconstruct  an  arbitrary  set  of  planes  in  the  breast  from  a  limited- 
angle  series  of  projection  images  when  the  x-ray  source  moves  along  an  arc  above  the  breast.  As  of  yet,  several 
tomosynthesis  algorithms  have  been  proposed  I3’4’10-22^  including  Shift-And-Add  (SAA)  [3,4],  Niklason  and 
colleagues’  publication  in  1997  of  a  tomosynthesis  method  with  the  x-ray  tube  moved  in  an  arc  above  the  stationary 
breast  and  detector  l0\  Wu  et  al.’s  report  in  2003  of  the  maximum  likelihood  iterative  algorithm  (MLEM)  to 
reconstruct  the  three-dimensional  distribution  of  x-ray  attenuation  in  the  breast  1 IJ,  and  the  filtered  back  projection 
(FBP)  algorithms  [12’13’i8"24l,  tuned-aperture  computed  tomography  (TACT)  reconstruction  methods  developed  by 
Webber  and  investigated  by  Suryanarayanan  et  a/[16],  algebraic  reconstruction  techniques  (ART)  F4-23-24^  filtered 
back  projection  (FBP)  '-4’13’18"21',  matrix  inversion  tomosynthesis  (MITS)  ^ 15,261  and  Gaussian  Frequency  Blending 
(GFB)  algorithm  to  combine  the  MITS  and  FBP  for  better  reconstruction  [22l 
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The  selection  of  optimal  acquisition  parameters  and  reconstruction  algorithm  plays  an  important  role  in  producing 
better  performance.  However,  there  are  several  factors  involved  in  this  task  and  some  of  them  are  not  individually 
independent.  An  effective  methodology  will  enable  one  to  optimize  selection  of  acquisition  parameters  and  compare 
algorithms  for  various  digital  breast  tomosynthesis  methods.  Recently,  Godfrey  et  al.  proposed  methods  of  MTF  and 
NPS  measurement  to  quantitatively  characterize  the  optimal  acquisition  parameters  selection  for  chest 
tomosynthesis  application  [6'7l  In  2005  and  2006,  we  also  published  the  impulse  response  and  NPS  analysis  for 
different  acquisition  parameters  and  reconstruction  algorithms  [3'4]. 

In  this  paper,  a  methodology  of  NEQ  (f)  analysis  is  proposed  to  optimize  and  compare  the  efficacy  of  tomosynthesis 
algorithms  and  imaging  acquisition  techniques  for  digital  breast  tomosynthesis.  It  combines  the  modulation  transfer 
function  (MTF)  of  system  signal  performance  and  the  noise  power  spectrum  (NPS)  of  noise  characteristics.  It 
enables  one  to  evaluate  the  performance  of  different  acquisition  parameters  and  algorithms  for  comparison  and 
optimization  purposes. 


2.  METHODS 


MTF2(f ) 

The  NEQ  (f)  is  defined  as:  NEO(  / )  = - - — .  The  modulation  transfer  function  (MTF)  and  the  normalized 

NNPS(f) 

noise  power  spectrum  (NNPS)  are  included  to  evaluate  the  performance  of  the  reconstruction  algorithms  and  the 

acquisition  techniques.  The  normalized  NPS  (NNPS)  is  defined  as:  NNPS(u,  v)  = - tomo  ^  ^  due  to  the 

gain 

logarithmic  transform  in  digital  tomosynthesis.  Hence,  the  NEQ  (f)  can  be  calculated  by: 


NEQ(f)  = 


gain2  -MTF2 {f) 

NPStomo(f ) 


.  In  digital  breast  tomosynthesis  reconstruction,  the  step  of  logarithmic  transform 


involved,  and  therefore,  the  determination  of  NEQ  (f)  is  more  similar  to  that  of  a  screen-film  system  than  to  a  linear 
digital  detector. 

2.1  MTF  measurement 

The  MTF  measurement  includes  two  parts:  1)  the  system  MTF  of  the  detector;  2)  the  reconstruction  MTF  associated 
with  specific  reconstruction  algorithm  and  acquisition  parameters. 


2.1.1  System  MTF  Measurement 

The  system  MTF  describes  the  measured  MTF  of  the  detector.  An  edge  method  was  applied  at  a  range  of  tube 
angles  to  see  if  there  is  any  difference  with  angle.  In  this  paper,  five  different  angular  positions  of  0fl,  ±15°,  and  ±25° 
of  the  X-ray  tube  were  selected  as  the  representative  angles. 

Table  1  shows  the  MW2  [SI  technique  we  used  for  the  system  MTF  measurement. 


Technique 

Acquisition 

Mode 

kV 

Spectrum 

Number  of 
views  (N) 

Total  mAs 

mAs  for  a  single 
projection 

MW2 

B0XD11  (slow 
speed) 

28 

W/Rh 

11 

303 

-21.5 

Table  1 .  System  MTF  measurement  technique 


53 


Figure  1  shows  the  set  up  of  the  system  MTF  measurement  experiment.  A  0.1mm  Pt-Ir  edge  was  placed  in  contact 
with  the  detector  and  oriented  at  a  1°~  3°  angle  with  respect  to  the  pixel  array.  Edge  images  were  then  acquired  with 
MW2  technique  and  2mm  Aluminum  filtration.  A  previously  published  routine  [8  9]  was  used  to  analyze  the  edge 
images  in  a  region  around  the  edge  (ROI=256*256  pixels,  pixel  size=85pm)  to  compute  the  presampled  system 
MTF. 


Edge  Device 

_ 

1°~3°  rotation 

_ 

Detector 


Chest  Wall 


Figure  1.  Set  up  of  system  MTF  measurements 

2.1.2  Reconstruction  MTF  measurement 

The  reconstruction  MTF  describes  the  calculated  MTF  associated  with  specific  algorithm  and  acquisition 
parameters.  A  dataset  of  tomosynthesis  projection  images  of  a  simulated  delta  function  with  the  specific  acquisition 
parameters  was  served  as  the  input  signal.  A  ray-tracing  simulation  method  was  used  to  project  the  single  delta 
function  onto  the  detector  to  simulate  the  tomosynthesis  sequence  of  projection  images.  Figure  2  shows  the  ray¬ 
tracing  method  to  simulate  the  projection  images  of  a  delta  function. 


Figure  2.  Ray-tracing  simulation 


In  Figure  2,  a  single  delta  function  at  z  distance  above  the  detector  was  simulated.  The  delta  function  is  projected 
onto  the  detector  at  P,  when  the  x-ray  tube  moves  to  T,  location.  Partial  pixel  sharing  was  performed  if  P,  falls  into 
non-integer  pixel  location.  With  this  ray-tracing  method,  for  a  specific  group  of  acquisition  parameters,  a  dataset  of 
projection  images  from  different  x-ray  tube  locations  can  be  simulated. 
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Projection  images  with  1/r  background  and  without  background  were  considered,  where  r  is  the  distance  between 
the  x-ray  source  7j  and  the  pixel  location  P,  on  each  projection  image.  Datasets  of  simulated  projection  images  with 
1/r2  background  and  without  background  were  separately  reconstructed  by  the  specific  reconstruction  algorithm.  The 
reconstruction  images  without  background  were  then  subtracted  from  the  reconstruction  images  with  background  to 
get  rid  of  effects  associated  with  the  1/r2  background  ray-tracing  simulation. 

The  reconstruction  MTF  can  be  computed  as  the  Fourier  Transform  of  the  point  spread  function  (PSF)  on  the 
reconstruction  plane  where  the  simulated  delta  function  is  located  (z  distance  above  the  detector). 

2.1.3  Total  MTF 

The  total  MTF  involves  two  parts  as  described  in  2.1.1  and  2.1.2.  It  can  be  calculated  as  the  multiplication  of 
measured  system  MTF  (2.1.1)  and  the  reconstruction  MTF  associated  with  reconstruction  algorithm  and  acquisition 
parameters  (2.1.2).  Linear  interpolation  was  performed  to  match  the  frequency  axes  of  the  system  MTF  and 
reconstruction  MTF  for  multiplication. 


2.2  NPS  measurement 

In  order  to  compute  the  noise  power  spectrum  (NPS)  for  specific  acquisition  parameters,  two  phantom  slabs  of 
BR12  (47%  water/  53%  adipose  equivalent)  for  a  total  of  4  cm  thickness  were  put  directly  on  the  detector  cover  to 
mimic  the  breast  tissue  equivalent  attenuation  and  scattered  radiation.  Ten  identical  tomosynthesis  sequences  of  flat 
images  with  the  phantom  slabs  on  the  detector  were  acquired.  The  tomosynthesis  sequences  were  then  reconstructed. 

Mean-subtracted  reconstruction  image  data  were  analyzed  and  compared  on  a  reconstruction  plane  with  the  same 
height  (z  distance  above  the  detector)  as  described  in  the  reconstruction  MTF  measurement  (2.1.2).  The  purpose  of 
studying  mean-subtracted  NPS  is  to  remove  fixed  pattern  noise,  including  structured  noise  and  system  artifacts.  A 
previously  published  method  [1 2]  was  applied  using  64  ROIs  of  size  128x128  to  examine  the  noise  response. 

2.3  Gain  factor  calculation 

In  order  to  normalize  NPS  and  fairly  compare  different  reconstruction  algorithms  and  acquisition  parameters,  the 
gain  factor  should  be  included  in  the  NEQ  (f)  calculation. 

For  an  x-ray  screen-film  imaging  system,  the  gain  factor  can  be  considered  as  gain  =  log10  (e)  ■  y  ,  where  y  is  the 

point  slope  of  the  film  density-log  x-ray  fluence  function  [28'29].  Figure  3  shows  a  typical  screen-film  nonlinear 
response  to  exposure.  Q  represents  the  number  of  incident  quanta  /  mm2.  From  this  nonlinear  response  curve,  one 
can  calculate  the  point  slope  y  for  a  specific  log(Q)  on  this  response  curve. 


With  digital  tomosynthesis,  the  gain  factor  varies  with  different  reconstruction  algorithms  and  acquisition 
techniques.  In  this  paper,  we  focus  on  the  relative  NEQ  (f)  methodology  of  a  single  algorithm  and  will  not  include 
the  gain  factor  calculation.  The  gain  factor  will  be  computed  in  future  work  where  we  will  actually  compare 
different  algorithms  against  one  another. 


Optical  Density 


Log  Relative  Exposure 


Figure  3.  Screen-film  characteristic  curve 
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2.4  NEQ  (f)  calculation 

gain2  ■  MTF2 (/) 

The  NEQ  (f)  can  be  calculated  by:  NEQ(f)  = - — — —  .  After  measuring  the  MTF,  NPS  and  gain 

NPStomo  (/) 

factor,  one  can  combine  them  together  to  get  the  NEQ  (f)  result.  Linear  interpolation  was  performed  to  match  the 
frequency  axes  of  the  NPS  and  MTF  curves.  In  this  paper,  the  relative  NEQ  (f)  along  two  directions  were 
investigated  separately,  including  the  x-ray  tube’s  movement  direction  and  the  direction  orthogonal  to  tube’s  motion 
direction. 


3.  RESULTS 

Here  we  give  an  example  of  NEQ  (f)  measurement.  A  prototype  Siemens  system  with  85  pm  pixel  size  was  used.  A 
point-by-point  Back  Projection  (BP)  reconstruction  algorithm  27  was  investigated.  Compared  with  traditional  shift- 
and-add,  point-by-point  BP  algorithm  demonstrates  improved  rendition  of  microcalcifications  in  the  direction 
perpendicular  to  the  tube  motion  direction,  thus  provides  sharper  appearance  of  calcifications  with  human  subject 
images  [27]. 

A  data  set  of  tomosynthesis  projection  images  of  a  delta  function  at  40  mm  above  the  detector  surface  plate  and  40 
mm  away  the  chest  wall  was  simulated  with  acquisition  parameters  of  25  projection  images  and  25  degrees  of  total 
angular  range  of  tube  movement.  The  simulated  tomosynthesis  sequence  was  reconstructed  by  point-by-point  BP. 
Figure  4  shows  the  impulse  and  MTF  results  on  the  reconstruction  plane  of  the  delta  function’s  location  (40  mm 
above  the  detector).  In  figure  4(a),  x  and  y  axes  represents  the  pixel  location  on  the  image.  In  figure  4(b),  the  v  axis 
represents  the  spatial  frequency  conjugate  to  the  direction  of  tube’s  motion  direction. 


Figure  4.  Point-by-point  BP  with  25  projections  and  ±15°  total  angular  range:  (a)  impulse  response;  (b) 


We  find  that  the  point-by-point  BP  is  an  effective  reconstruction  method  for  digital  breast  tomosynthesis.  The  MTF 
curve  is  smooth  and  the  impulse  response  appears  to  be  sharp  at  the  reconstruction  plane  where  the  simulated  delta 
function  is  located  (figure  4). 


Figure  5  shows  the  edge-method  measured  system  MTF  when  the  edge  center  was  about  4cm  away  from  the  chest 
wall  and  the  X-ray  tube  was  at  angles  of  0°,  ±15°,  ±25°.  Three  datasets  were  measured  and  averaged  for  each  angular 
location  of  the  x-ray  tube.  The  MTF  varied  little  with  different  angles.  The  MTF  curve  of -25"  is  a  little  lower  than 
that  of  other  angles.  This  may  be  caused  by  x-ray  tube’s  motion  and  velocity  difference  at  the  -25°  location,  which  is 
the  beginning  position  of  the  x-ray  tube  during  the  tomosynthesis  sequence. 
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Figure  5.  System  MTF  at  different  angles  by  edge  method  measurement 


system  MTF 

0 _ _  -H-  total  MTF _ J. _ i _ _ 

0  1  2  3  4  5  6 

Frequency  (cycles/mm) 


total  MTF 

O' - 1 - 1 - * - 1 - 1 - 1 

0  1  2  3  4  5  6 

cycles/mm 


(a) 


(b) 


Figure  6.  Reconstruction,  system,  and  total  MTFs:  (a)  direction  orthogonal  to  tube’s  motion  direction  (b)  tube’s  motion  direction 


Figure  6  shows  the  averaged  system  MTF,  reconstruction  MTF  (point-by -point  BP  algorithm,  25  projection  images 
and  ±12. 5n  total  tube  angular  movement),  and  the  total  MTF  as  the  multiplication  of  above  two  MTFs.  The 
reconstruction  MTF  is  much  higher  than  the  system  MTF. 

For  measurement  of  NPS,  a  W/Rh  spectrum  and  28  KV  were  used  with  a  cumulative  tube  output  of  358  mAs. 
Mean-subtracted  NPS  of  reconstructed  planes  at  the  same  location  of  40mm  above  the  detector  surface  plate  was 
computed  as  shown  in  figure  7.  Figure  8  shows  the  relative  NEQ  (f)  results  without  gain  factor  calculation.  Further 
theoretical  and  experimental  studies  will  be  done  to  investigate  the  gain  factor  for  each  reconstruction  algorithm  and 


57 


acquisition  technique.  Results  along  both  tube’s  movement  direction  and  direction  orthogonal  to  the  tube’s  motion 
are  shown  in  figures  7  and  8.  The  NPS  and  relative  NEQ  (f)  performances  along  both  directions  were  similar. 


cycles/mm  cycles/mm 

Figure  7.  Mean- subtracted  NPS  Figure  8.  Relative  NEQ  (f) 

4.  CONCLUSION 

The  NEQ  (f)  is  an  image  quality  metric  that  combines  both  signal  and  noise  properties  to  compare  different 
acquisition  techniques  and  reconstruction  algorithms  in  tomosynthesis.  In  this  paper,  MTF  and  NPS  were  evaluated 
by  experiments,  simulation,  and  the  application  of  several  published  routines.  This  provides  empirical  results  for 
comparing  and  selecting  optimal  tomosynthesis  acquisition  parameters  and  reconstruction  algorithms,  and  enables 
elucidation  of  imaging  geometry  factors  in  directions  parallel  and  orthogonal  to  tube  motion. 
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